; docformat = 'rst'
;
;+
;
; :Purpose:
;   Regrid GAUGE flux data for use in MOZART
;
; :Inputs:
;
; :Requires:
;
; :Keywords:
;
; :Outputs:
;
; :Example::
;
; :History:
; 	Written by: Matt Rigby, University of Bristol, Nov 25, 2013
;
;-
pro gauge_n2o_flux_mozart_regrid

  ;Get flux data
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  species='N2O'
  Title='N2O_emissions' ;Global "title" attribute for netCDF file
  Author='Matt Rigby'

  ncName=XXXX ;output file name

  ;Input file name
  YOU will want to loop this first bit over all the input fluxes and add them together
  
  flux_fname='/data/GAUGE/N2O_flux/N2O_flux_natsoil_global_1990_2008.nc'

  lon_ncname='lon'
  lat_ncname='lat'
  flux_ncname='flux'

  q=mr_ncget(flux_ncname, filename=flux_fname) ;kg/m2/s
  q_lat=mr_ncget(lat_ncname, filename=flux_fname)
  q_lon=mr_ncget(lon_ncname, filename=flux_fname)
  date=mr_ncget('date', filename=flux_fname)
  datesec=mr_ncget('datesec', filename=flux_fname)

  q_area=mr_areagrid(q_lon, q_lat)

  mr_grid_orient, q_lon, q_lat, q
  
  time=mr_mozdate(date)
  time+=double(datesec)/24.d/3600.d

  
  ;Get MOZART grid
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  
  MOZART_fname='/data_rdsf/MOZART/MERRA/MERRA_19x2_20000101.nc'
  lon=mr_ncget('lon', filename=MOZART_fname)
  lat=mr_ncget('lat', filename=MOZART_fname)
  area=mr_areagrid(lon, lat)
  
  
  ;Define output array
  q_regrid=fltarr(n_elements(lon), n_elements(lat), n_elements(time))  

  ;Martin Shultz's regridding routine.
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  
  ;Define old grid
  nlon = n_elements(q_lon)
  incr = (max(q_lon) - min(q_lon))/Double(nlon-1)
  ms_lon = Obj_New('MGS_LonVector', $
    start=min(q_lon),incr=incr,nval=nlon, $
    type='regular')
  nlat = n_elements(q_lat)
  incr = (max(q_lat) - min(q_lat))/Double(nlat-1)
  ms_lat = Obj_New('MGS_LatVector', $
    start=min(q_lat),incr=incr,nval=nlat, $
    type='regular')
  o = Obj_New('MGS_RGrid', name='Source Grid', $
    Longitude=ms_lon, Latitude=ms_lat)
  if total(o->getlon() - q_lon) gt 0.1 then stop
  if total(o->getlat() - q_lat) gt 0.1 then stop
  
  ;Define new grid
  LonSize = n_elements(lon)
  lon_incr = (max(lon) - min(lon))/Double(LonSize-1)
  ms_lon = Obj_New('MGS_LonVector', $
    start=min(lon),incr=lon_incr,nval=LonSize, $
    type='regular')
  LatSize = n_elements(lat)
  lat_incr = (max(lat) - min(lat))/Double(LatSize-1)
  ms_lat = Obj_New('MGS_LatVector', $
    start=min(lat),incr=lat_incr,nval=LatSize, $
    type='regular')
  n = Obj_New('MGS_RGrid', name='Target Grid', $
    Longitude=ms_lon, Latitude=ms_lat)
    
  if total(n->getlon() - lon) gt 0.1 then stop
  if total(n->getlat() - lat) gt 0.1 then stop
  
  for ti=0, n_elements(time)-1 do begin
    ;Regrid
    q_regrid[*, *, ti]=regrid(q[*, *, ti], oldgrid=o, newgrid=n, /use_grid, /area_weighted) ;kg/m2/s
    print, string(float(ti)/float(n_elements(time))*100., format='(I)'), '% complete. ' + $
      string(round(total(q_regrid[*, *, ti]*q_area)/total(q[*, *, ti]*area)-1.)*100., format='(I)'), '% error'
  endfor
  
  ;Clean
  Obj_Destroy, o
  Obj_Destroy, n


  CONVERT TO MOLECULES/CM2/S

  ;Now write in MOZART netCDF file
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
  fID=ncdf_create(ncName, /clobber)
  
  LonDID  = NCDF_DIMDEF(fID, 'lon', long(n_elements(lon)))
  LatDID  = NCDF_DIMDEF(fID, 'lat', long(n_elements(lat)))
  TimeDID  = NCDF_DIMDEF(fID, 'time', /unlimited)
  
  ; Create a variable to hold the data
  lonID = NCDF_VARDEF(fId,'lon', [LonDID], /DOUBLE)
  latID = NCDF_VARDEF(fId,'lat', [LatDID], /DOUBLE)
  dateID = NCDF_VARDEF(fId,'date', [TimeDID], /LONG)
  timeID = NCDF_VARDEF(fId,'time', [TimeDID], /DOUBLE)
  qID = NCDF_VARDEF(fId, strcompress(strlowcase(species) + 'emissions', /remove_all), [LonDID, LatDID, TimeDID], /DOUBLE)
  
  ;Create some global attributes
  NCDF_ATTPUT, fid, /Global, "Title", title
  NCDF_ATTPUT, fid, /Global, "Author", author
  NCDF_ATTPUT, fid, /Global, "Created", systime()
  
  ; Create some attributes (about our variables)
  NCDF_ATTPUT, fId, lonId, "long_name", "longitude"
  NCDF_ATTPUT, fId, lonId, "units", "degrees_east"
  
  NCDF_ATTPUT, fId, latId, "long_name", "latitude"
  NCDF_ATTPUT, fId, latId, "units", "degrees_north"
  
  NCDF_ATTPUT, fId, qId, "long_name", Strcompress(strupcase(species) + " emission rate")
  NCDF_ATTPUT, fId, qId, "units", "molecules/cm2/s"
  
  NCDF_ATTPUT, fId, TimeID, "long_name", "Time"
  NCDF_ATTPUT, fId, TimeID, "units", "days since 0000-01-01 00:00:00"
  
  NCDF_ATTPUT, fId, DateID, "long_name", "Date"
  NCDF_ATTPUT, fId, DateID, "units", "YYYYMMDD"
  
  ; Leave definition mode and enter data write mode
  NCDF_CONTROL, fId, /ENDEF
  
  ; Write the data
  NCDF_VARPUT, fId, lonId, lon
  NCDF_VARPUT, fId, latId, lat

  NCDF_VARPUT, fId, dateID, mr_mozdate(time)
  NCDF_VARPUT, fId, timeID, time - julday(1, 121, 1, 0)
  NCDF_VARPUT, fId, qId, q
  
  ; Done
  NCDF_CLOSE, fId


  ;FOR 3D fields, I'd use mr_hybrid_coords() to interpret the hyam and hybm parameters
  ; Then use regrid() to regrid based on pressure and lon/lat...

end